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ABSTRACT 

The turbulent stresses that lead to angular momentum transport in accretion discs 
have often been treated as resulting from an isotropic effective viscosity, related to 
the pressure through the alpha parametrization of Shakura & Sunyaev. This simple 
approach may be adequate for the simplest aspects of accretion disc theory, and was 
necessitated historically by an incomplete understanding of the origin of the turbu- 
lence. More recently, Balbus & Hawley have shown that the magnetorotational in- 
stability provides a robust mechanism of generating turbulent Reynolds and Maxwell 
stresses in sufficiently ionized discs. The alpha viscosity model fails to describe numer- 
ous aspects of this process. This paper introduces a new analytical model that aims 
to represent more faithfully the dynamics of magnetorotational turbulent stresses and 
bridge the gap between analytical studies and numerical simulations. Covariant evo- 
lutionary equations for the mean Reynolds and Maxwell tensors are presented, which 
correctly include the linear interaction with the mean flow. Non-linear and dissipative 
effects, in the absence of an imposed magnetic flux and in the limit of large Reynolds 
number and magnetic Reynolds number, are modelled through five non-linear terms 
that represent known physical processes and are strongly constrained by symmetry 
properties and dimensional considerations. The resulting model explains the develop- 
ment of statistically steady, anisotropic turbulent stresses in the shearing sheet, a local 
representation of a differentially rotating disc, in agreement with numerical simula- 
tions. It also predicts that purely hydrodynamic turbulence is not sustained in a flow 
that adequately satisfies Rayleigh's stability criterion. The model is usually formally 
hyperbolic and therefore 'causal', and guarantees the realizability of the stress tensors. 
It should be particularly useful in understanding the dynamics of warped, eccentric 
and tidally distorted discs, non-Keplerian accretion flows close to black holes, and a 
variety of time-dependent accretion phenomena. 
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1 INTRODUCTION 

The magnetorotational instability is one of the most im- 
portant instabilities in astrophysical fluid dynamics. It ap- 
plies to a differentially rotating, electrically conducting fluid 
in which the angular velocity decreases in magnitude away 
from the axis of rotation. 1 In the presence of a weak mag- 
netic field of arbitrary configuration, such a flow is sub- 
ject to a dynamical instability, with a growth rate compa- 
rable to the shear rate of the flow (Velikhov 1959; Chan- 
drasekhar 1960; Fricke 1969; Acheson 1978; Balbus & Haw- 
ley 1991, 1992; Papaloizou & Szuszkiewicz 1992; Balbus 



1 In the absence of entropy gradients. Otherwise, modified 
versions of the Hoiland stability criteria exist (Papaloizou & 
Szuszkiewicz 1992; Balbus 1995). 



1995; Foglizzo & Tagger 1995; Ogilvie & Pringle 1996; 
Terquem & Papaloizou 1996). 

The principal application of the magnetorotational in- 
stability is to accretion discs, in which the profile of angular 
velocity is fixed by Kepler's third law. The non-linear de- 
velopment of the instability leads to sustained magnetohy- 
drodynamic (MHD) turbulence, which transports angular 
momentum outwards in a vain attempt to neutralize the 
destabilizing gradient of angular velocity (Hawley, Gammie 
& Balbus 1995; Brandenburg et al. 1995; Stone et al. 1996; 
Balbus & Hawley 1998). Owing to the generality of the con- 
ditions for instability, however, further applications exist to 
stellar interiors and other astrophysical objects, and possi- 
bly to geophysical situations. Laboratory experiments are 
also in progress (Ji, Goodman & Kageyama 2001). 

In one sense, the magnetorotational instability is gen- 
erally considered to have 'solved' the problem of the effec- 
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tive viscosity of accretion discs, through providing a robust 
means of angular momentum transport. 2 However, the re- 
alization of the importance and ubiquity of the instability 
has, in some ways, made matters more difficult. A strict view 
might be taken that any calculation of an accretion disc must 
now involve a numerical simulation of three-dimensional 
MHD turbulence. These simulations are demanding to ex- 
ecute, complicated to analyse and still restricted in their 
scope and physical content. 

To illustrate the limitations of the direct numerical ap- 
proach, let us consider how the computational requirements 
scale with the parameters of the problem being studied. 
Most simulations conducted so far have been based on the 
shearing box (Hawley et al. 1995), which is a local represen- 
tation of a differentially rotating disc. The simulated volume 
is comparable to H 3 , where H is the vertical scale-height or 
semi-thickness of the disc, and the models can be run for 
many times the dynamical (orbital) time-scale. Let us ac- 
cept, for the purposes of this argument, that recent versions 
of such simulations (e.g. Miller & Stone 2000) have achieved 
an adequate fidelity to physical reality (although questions 
indeed remain concerning numerical resolution and conver- 
gence, vertical boundary conditions, thermal and radiative 
physics, etc.). How much more difficult would it be to simu- 
late a global model of a realistic thin disc over many times 
the viscous time-scale, as would be required in a typical 
application? Suppose that the disc has a constant angu- 
lar semi-thickness H/r and is simulated in a conical wedge 
in spherical polar coordinates (r, 0, (f>) with logarithmic grid 
spacing in r between n n and r out , and uniform spacing in 
9 and <j>. To achieve a resolution comparable to that of a 
shearing box, the number of grid zones must increase by a 
factor of order (r/H) ln(r ou t/rin) in r, and a factor of or- 
der r/H in <j>. The global viscous time-scale exceeds the 
dynamical time-scale in the inner part of the disc (where 
the Courant condition on the time-step is most severe) by 
a factor of order (l/a)(r/-ff) 2 (r ou t/n n ) 3 ' /2 , where a is the 
viscosity parameter of Shakura & Sunyaev (1973). There- 
fore the computational requirements increase by a factor of 
order 

Reasonable estimates of this factor for discs around young 
stars, in cataclysmic variables, and in X-ray binaries, respec- 
tively, are 10 13 , 10 11 and 10 15 , or even larger. (These num- 
bers are reduced if the inner part of the disc is truncated 
by a stellar magnetosphere.) For this reason, although some 
aspects of large-scale disc dynamics can now be addressed 
in direct numerical simulations (e.g. Hawley 2001), global 
studies of realistic thin discs are currently inaccessible. 

How then are we to study situations such as the 
thermal-viscous instability of discs in cataclysmic variables, 
the dynamics of warped or eccentric discs, tidally dis- 
torted discs in binary stars and protoplanetary systems, 
non-Keplerian accretion flows close to black holes, or the 

2 Nevertheless, considerable interest and debate remain regarding 
the operation of the magnctorotational instability, or alternative 
mechanisms of angular momentum transport, in weakly ionized 
discs, e.g. around young stars and in the quiescent state of dwarf 
novae. 



propagation of waves in discs? The traditional approach, 
of course, has been to treat the turbulence as an effective 
viscosity given by the alpha parametrization of Shakura & 
Sunyaev (1973). Despite suspicions that this treatment is 
likely to be inadequate in various respects, very few studies 
have aimed specifically at testing this issue (Abramowicz, 
Brandenburg & Lasota 1996; Torkelsson et al. 2000) and 
it is still unclear in what respects magnetorotational tur- 
bulence does or does not behave as an effective viscosity. 
In a standard accretion disc the rate-of-strain tensor has a 
dominant r^i-component and an isotropic viscous model pre- 
dicts that the rcS-component of the turbulent stress would 
similarly dominate. This is not true of magnetorotational 
turbulence, in which other stress components (e.g. cf>(j>) are 
larger. However, for the simplest problems such as the evolu- 
tion of the surface density of a standard accretion disc, these 
additional stress components play essentially no role in the 
dynamics. The alpha model is probably perfectly adequate 
in these circumstances, where one dimensionless number, a, 
is parametrizing a single physical quantity, the vertically in- 
tegrated stress coefficient T r ^. 3 

In situations such as warped, eccentric, tidally dis- 
torted and non-Keplerian discs there is a large-scale, pos- 
sibly slowly varying velocity field that is different from that 
in a standard disc. From the perspective of a local observer, 
however, the flow resembles a shearing-box model, but with 
a non-standard (and possibly slowly varying) velocity gradi- 
ent tensor. Stress components other than the ^-component 
can play a significant role in the dynamics of these systems, 
but whether the full turbulent stress tensor resembles a vis- 
cous stress in such situations is not understood. In princi- 
ple, shearing-box studies could be used to formulate a more 
faithful analytical or semi-analytical model of the relevant 
properties of the turbulence, which could then be used in 
studies of the large-scale behaviour of accretion discs. This 
would bridge the gap between the local dynamics (on scales 
from the dissipation length up to H) and the large-scale dy- 
namics, without requiring colossal advances in computation. 
The formation and analysis of such a model is also likely to 
result in a better physical understanding. 

In view of the formidable complexity of turbulent flows, 
such a programme might appear overarnbitious. A com- 
plete theory of turbulence does not exist, because turbulent 
flows have an unlimited number of statistical properties that 
cannot be calculated from first principles. However, in the 
present context, our attention is restricted to the question 
of how the mean turbulent stress in a patch of the disc re- 
sponds to changes in the large-scale velocity gradient. As will 
be made clear, certain linear aspects of this problem can be 
deduced directly from the MHD equations; other non-linear 
aspects cannot be treated rigorously and require a closure 
model that is, nevertheless, strongly constrained by symme- 
try properties and dimensional considerations. Although the 
model affords only an imperfect description of turbulence, its 
physically motivated nature and its connection to the MHD 
equations ensure that it performs more faithfully than the 



' 5 It may not be adequate for predicting the vertical distribution 
of energy dissipation, and therefore the vertical structure of the 
disc. 
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alpha viscosity model, and this is borne out by numerous 
test problems. 

The remainder of this paper is organized as follows. In 
Section 2 I introduce a conceptual model system for magne- 
torotational turbulence, and make some simple arguments 
concerning the saturation of the turbulence. I develop the 
equations for the mean Reynolds and Maxwell tensors as far 
as possible analytically in Section 3, and then discuss the re- 
quirements of a non-linear closure model. A specific model 
satisfying these principles is presented and explained physi- 
cally in Section 4. The outcome of the model in the shearing 
sheet, under a variety of conditions, is investigated in Sec- 
tion 5. In Section 6 the effects of a compressible mean flow 
are incorporated and the total energy budget is explained. 
Some preliminary applications of the model are worked out 
in Section 7. In Section 8 a comparison is made with ex- 
isting models in engineering and in astrophysics, and also 
with published numerical simulations. A summary follows 
in Section 9. 



2 A MINIMAL MODEL SYSTEM FOR 

MAGNETOROTATIONAL TURBULENCE 

The magnetorotational instability has certain minimal re- 
quirements: rotation and shear (in the correct relative ori- 
entation), electrical conductivity and a seed magnetic field. 
Let us consider the simplest conceptual model system for 
studying magnetorotational turbulence: the incompressible 
shearing sheet (Goldreich & Lynden-Bell 1965). Let (x,y,z) 
be Cartesian coordinates in a frame of reference rotating 
with uniform angular velocity $7 = fi e z . An incompressible 
fluid of uniform density p, kinematic viscosity v and mag- 
netic diffusivity r) has a uniformly shearing velocity field 
characterized by a uniform velocity gradient tensor Vu. In 
a standard disc this is of the form u = —2Axe y , where 
A is Oort's first constant. The x, y and z directions corre- 
spond to the radial, azimuthal and vertical directions from 
the perspective of a corotating observer in a thin, differen- 
tially rotating disc. The flow is unbounded in the horizontal 
(xy) plane but of bounded vertical extent, with either 'phys- 
ical' or periodic boundary conditions. The distance between 
the vertical boundaries defines a characteristic length-scale 
L = L z which, in a real disc, is related to the vertical scale- 
height or thickness. 

As defined above, the system involves just three dimen- 
sionless parameters: the Rossby number Ro = A/Q,, the 
Reynolds number Re = L 2 Q,/v and the magnetic Prandtl 
number Pm = v jr\. In a Keplerian disc, the Rossby number 
is fixed at Ro = 3/4. The microscopic diffusivities v and r\ 
are included in this conceptual model only to regularize the 
system by allowing for dissipation and irreversibility. The 
Reynolds number may be considered to be arbitrarily large. 

Apart from the assumption of incompressibility and the 
inclusion of microscopic diffusivities, this system is equiva- 
lent to the shearing box studied by Hawley et al. (1995) in 
the limit that the horizontal scales L x and L y of the box 
tend to infinity while retaining a finite vertical scale. Note 
that the limit L x ,L y S> L z corresponds to the physical sit- 
uation in a thin disc. It is natural to make the following 
assumptions, which are not contradicted by existing numer- 
ical simulations: 



• the macroscopic statistical properties of the turbulence 
do not depend on the details of the dissipative mechanisms, 
and are therefore independent of Re and Pm in the limit 

Re — > oo; 

• the macroscopic statistical properties of the turbulence 
are bounded and well defined in the limit L x ,L y — > oo. 

The system is horizontally homogeneous, in the sense 
that any point in the xy-plane is equivalent (modulo a 
Galilean boost) . If periodic vertical boundary conditions are 
imposed, the system is also vertically homogeneous. It is 
then possible for the system to develop statistically steady 
and homogeneous (although anisotropic) turbulence. 

If magnetorotational turbulence is initiated in this sys- 
tem, the properties of the saturated state are strongly con- 
strained by dimensional considerations. The RMS turbulent 
velocity fluctuation, (u 2 ) 1 ^ 2 , for example, must equal LO 
times a dimensionless coefficient of order unity, if the as- 
sumed independence of Re and Pm is correct. Therefore the 
vertical scale of the system plays an essential role in the non- 
linear saturation of the turbulence, presumably by limiting 
the size of coherent structures ('eddies'). In an unbounded 
system with L — * oo, the turbulence would presumably grow 
indefinitely without saturation. It is sometimes suggested 
that turbulent motions and magnetic fields in an accretion 
disc are limited by shock formation or magnetic buoyancy, 
but there is little evidence of this in numerical simulations 
of magnetorotational turbulence, and neither limiting mech- 
anism operates in the incompressible shearing sheet. 4 



3 STRESS EQUATIONS AND CLOSURES 
3.1 Basic equations 

Consider an incompressible fluid of uniform density p, kine- 
matic viscosity v and magnetic diffusivity r\. In the Carte- 
sian tensor notation, the equation of motion in a frame of 
reference rotating with uniform angular velocity fi; is 

(d t + Ujdj)ui + 2eij k 0.jU k = -<9iII + bjdjbi + vdjjUi, (2) 

where 

bi = { lM)P y 1/2 B l (3) 
is the Alfven velocity and 

n = £ + \b 2 + $ (4) 
p 

is the modified pressure including the gas pressure, the mag- 
netic pressure and the gravitational and centrifugal poten- 
tials. The induction equation is 

(d t + Ujd 3 )bi = bjdjUi + rjd n b t , (5) 

and the velocity and magnetic fields satisfy the solenoidal 
conditions, 

diUi = dibi = 0. (6) 

Separate the varying quantities into mean and fluctu- 
ating parts, e.g. 

4 Numerical simulations of an incompressible shearing box us- 
ing a spectral method (unpublished work by the author) produce 
results similar to those of Hawley et al. (1995). 
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Ui =Ui + u'i, (7) 

where Ui = (ui) and (u[) = 0. The angle brackets denote a 
suitable averaging procedure such as a spatial, temporal or 
ensemble average. The mean forms of the above equations 
are then 



(d t + Ujdj)ui + 2e ijk QjU k = — diU + bjdjbi + vdjjUi 



+dj(Mij -Rij), 


(8) 


(dt + Ujdj)bi = bj9j"Ui + ridjjbi + djFij, 


(9) 


d^i = dibi = 0, 


(10) 


where 




Mi j = b'ib'j, 


(11) 


JXij — UiUj , 


(12) 




(13) 



are the Maxwell, Reynolds and Faraday tensors. The prob- 
lem at hand is how to determine the mean quantities Mij, 
Rij and Fij and thereby close the system of equations. 5 

Unless the mean velocity field acts as a dynamo, the 
mean magnetic field can be sustained only by the mean 
Faraday tensor, which represents the turbulent electromo- 
tive force or 'alpha effect'. For simplicity, it will be assumed 
in the following that there is no mean magnetic field and no 
large-scale dynamo, i.e. bi = Fij — 0. The remaining prob- 
lem is to determine the Maxwell and Reynolds stress tensors 
that appear in the mean equation of motion. 

In the case bi = Fj = 0, the induction equation reads 

(d t + Ujdj)b'i + u'jdjb'i = b'jdjUi + b'jdjul + r/djjbi- (14) 

From this, an exact equation for the mean Maxwell stress 
can be obtained, in the form 

(d t + Ukdk)Mij — M ik dkUj — MjkdkUi 

= (b'idkFjk + b'jdkFik) + ^(b'idkkb'j + b'jd kk b'i). (15) 

The left-hand side of this equation represents the linear dy- 
namics of the Maxwell tensor. It shows how the tensor is 
advected by the mean velocity field and 'stretched' by gra- 
dients in the mean velocity. The right-hand side contains 
difficult terms of two types: those involving triple correla- 
tions of fluctuating quantities, and those involving the mi- 
croscopic diffusion process. These are both 'non-linear' ef- 
fects; although Ohmic diffusion is a linear process, when 
the magnetic Reynolds number is large the diffusive terms 
can be significant only when a turbulent cascade has forced 
structure to appear on the dissipative scales. 

There is little hope of finding a rigorous closure scheme 
for this equation. Instead, the approach adopted here will 
be to write the equation in the form 

(d t + u k dk)Mij — MikdkUj — MjkdkUi = • ■ ■ , (16) 

retaining the exact form of the linear terms, and then inves- 
tigate simple closure models for the non-linear terms on the 
right-hand side. 

5 For the high-Reynolds number situations considered in this pa- 
per, the diffusive terms involving v and r\ in the mean equations 
may be neglected. 



A similar approach applied to the fluctuating part of 
the equation of motion leads to the exact equation 

(dt + Ukd k )Rij + RikdkUj + RjkdkUi 

+2e J ki^kRu + 2eikiSlkRji = -{u'idjU' + u'jdiU') 

+ {u'idk(Mjk - Rjk) +u' J dk(M i k - Rik)) 

+v(u' i dkku' :j + u'jdkku'i) (17) 

for the mean Reynolds tensor. The terms on the right-hand 
side require a closure model, but the form of the linear terms, 

(d t + u k d k )Rij + RikdkUj + RjkdkUi 

+2ejkiQkRu + 2etkiflkRji = • • • , (18) 

may be retained. Note that the Maxwell and Reynolds ten- 
sors interact with the mean velocity gradient in different 
ways, and only the Reynolds tensor is affected by rotation. 

3.2 The pressure— strain correlation 

In homogeneous turbulence, the term involving II' on the 
right-hand side of equation (17) has the alternative form 

(Yl'(dju'i + diu'j)), (19) 

and is known as the pressure-strain correlation. The iden- 
tification of suitable closures for this term, even in purely 
hydrodynamic turbulence, has been a matter of consider- 
able controversy (e.g. Speziale 1991) leading to highly elab- 
orate, but still imperfect, models (e.g. Sjogren & Johansson 
2000). An expression for II' can be obtained by taking the 
divergence of the fluctuating part of the equation of mo- 
tion. If the spectrum of the turbulence is known, part of the 
pressure-strain correlation can then be expressed in terms of 
Rij, while part is undoubtedly non-linear. However, as the 
spectrum itself is determined through non-linear dynamics, 
one may take the view that the entire pressure-strain corre- 
lation is a non-linear term for which a non-deductive closure 
must be proposed. 

In the following, the fluctuating quantities will not be 
referred to, and the bars on Ui, Mij and Rij will be omitted. 

3.3 Linear dynamics in the shearing sheet 

In the shearing sheet the angular velocity is = Q e z 
and the only non- vanishing component of the mean velocity 
gradient is d x u y — —2A. The linearized equations for the 
Maxwell and Reynolds stresses then have the form 

d t M xx = 0, 

d t M xy + 2AM XX = 0, 

d t M xz = 0, 

d t Myy + AAM X y = 0, 

d t M yz + 2AM XZ = 0, 
d t M zz = 0, 

dtR xx — 4:Q.R X y = 0, 

dtR xy + 2(f2 — A)R XX — 2Q.R yy = 0, 
dtR xz — 2flR yz = 0, 

dtRyy + 4(fi - A)R xy = 0, 

d t Ryz + 2(Q - A)R XZ = 0, 
d t R zz = 0. 
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The Maxwell and Reynolds tensors are decoupled. The gen- 
eral solution for the magnetic component is 

M xx = M xx0 , 

M xy = M xy0 - 2M xx0 At, 

M xz — M xz o, 

M yy = M yy0 - ±M xy0 At + AM xx0 A 2 t 2 , 
M yz = M yz0 - 2M xz0 At, 
M zz = M zz0 , 

and allows for algebraic growth through the shearing of mag- 
netic field lines. The solution for the kinetic component de- 
pends on the Rayleigh discriminant, or squared epicyclic fre- 
quency, 

K 2 =4n(n-A), (20) 

of the system. In the absence of rotation, the system is 
Rayleigh-neutral and the solution for Rij allows for alge- 
braic growth (it resembles the solution for Mij, but with 
the sign of A reversed). When k 2 > 0, the solution is oscil- 
latory, as can be seen by reducing the problem to equations 
of the form 

dttR xy +4k 2 R xv = 0, 
d u Rxz + 2kR xz = 0, 
duRzz = 0. 

So R xx , R xy and R yy oscillate at twice the epicyclic fre- 
quency, R xz and R yz at the epicyclic frequency, and R zz is 
constant. On the other hand, when k 2 < 0, Rayleigh's crite- 
rion for stability is violated and exponential growth occurs. 

3.4 Requirements of a non-linear closure model 

Whenever the linear dynamics indicates unbounded growth 
or undamped oscillations, the non-linear terms will deter- 
mine the eventual outcome. The approach adopted in this 
paper is to explore the consequences of simple closure mod- 
els of the form 

V {1) Rij = ^{Ri^Ma,...), 

V^Mij = ^{Ri^Ma,...), (21) 

where the operators T> are defined by the left-hand sides of 
equations (18) and (16) respectively, and the quantities Tij 
are non-linear tensorial functions of their arguments. The 
dots represent the parameters of the problem, on which the 
functions Tij may depend. 

Such a model ought to satisfy the following fundamental 
principles. 

• There should be no 'source terms'. Rij = Mij = 
should always be a solution, although it may be unstable. 

• An unmagnetized state, Mij = 0, should always be a 
solution even when Rij / 0, although it may be unstable to 
dynamo action. 

• The non-linear terms on the right-hand side should be 
manifestly covariant, as the left-hand sides are. 

• The positive semi-definite nature of the tensors Rij and 
Mij , implicit in their definition, should be preserved by the 
model. 6 This ensures not only that the turbulent kinetic 

6 A real symmetric n X n matrix (or second-rank tensor) Aij 



and magnetic energy densities remain non- negative, but also 
that the stress tensors can be realized by genuine velocity 
and magnetic fields. 

The following two additional requirements appear plausible, 
although they cannot be regarded as fundamental principles. 

• The non-linear terms should not refer to the angular 
velocity fl or the velocity gradient Vu, because their effects 
are fully represented in the linear terms. 

• The non-linear terms should not refer to the mi- 
croscopic diffusion coefficients, because one has assumed 
asymptotic independence of Re and Pm in the limit Re — > 
oo. 

The non-linear terms are also strongly constrained by di- 
mensional considerations. Both Rij and Mij, as defined in 
this incompressible system, have the dimensions of velocity- 
squared. A quantity with the dimensions of the rate of 
change of Rij cannot be formed from Rij and My alone. The 
only other physical quantity that can appear in the functions 
T , according to the above assumptions, is the length-scale 
L. In that case the non-linear terms must have the form 

J*? =L- 1 g%\R i j,M ij ), (22) 

where the quantities Qij are homogeneous functions of de- 
gree 3/2. 

4 A SIMPLE MODEL AND ITS PROPERTIES 

4.1 Statement of the model 

A simple model based on the above principles is as follows. 

(d t + u k d k )Rij + R ik d k uj + Rj k d k Ui 
+2ej k iQ k Ru + 2ei k iQ k Rji 
= -C 1 L- 1 R 1/2 R l , - C 2 L- 1 R 1/2 (R lj - ^RSij) 

+C 3 L- 1 M 1/2 M^ - C4,L~ 1 R~ 1 / 2 MRij, (23) 

(d t + u k d k )Mij — M ik d k Uj — M jk d k Ui 

= dL^R-^MRij - (C 3 + Cr,)L- 1 M 1/2 M lJ , (24) 

where R = R it and M = Ma are the traces of the Reynolds 
and Maxwell tensors, and Ci , . . . , Cs are positive dimen- 
sionless coefficients of a universal nature. 

In the following subsections the physical origin and im- 
plications of the various terms in the model will be made 
clear. 

4.2 Hydrodynamic decay and return to isotropy 

For freely decaying hydrodynamic turbulence in the absence 
of rotation and a mean flow, one might start with a model 
of the form 

is said to be positive semi-definite if AijXiXj > for all real 
n-component vectors X{, or, equivalently, if all the eigenval- 
ues of Aij are non-negative. This condition is relevant because 
RijX % Xj = (u'^XiXj) = {(u' ■ X) 2 ) > 0, and similarly for 
Mij. Note that off-diagonal components such as R xy can have 
either sign. 
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d t R l3 • = -C 1 L- 1 R 1/2 R lj . (25) 

The turbulent kinetic energy density ^R then decays mono- 
tonically and non-linearly according to 

d t R= -dL^R^ 2 . (26) 

The functional form of the right-hand side is dictated by 
dimensional considerations, as discussed above, and is also 
suggested by the form of the triple correlations that appear 
in the exact equation (17). The predicted decay of energy, 
R oc t~ 2 as t — ► oo, is faster than the behaviour R oc t~ x 
observed in turbulence generated by the passage of a stream 
of fluid through a grid (Batchelor 1953). The reason for this 
is that the length-scale of the energy-containing eddies is 
constrained by the constant L in the present model, while 
it expands continuously during the decay of grid-generated 
turbulence. 

It is well known that freely decaying, anisotropic hy- 
drodynamic turbulence also exhibits a tendency to return 
to isotropy. This suggests the addition of a traceless term 
that redistributes energy among the components of Rij, i.e. 

dtRij = ~C\L- 1 R 1/2 R lJ - C2L- 1 R 1/2 {R ij - ±R5ij). (27) 

In this model, if initially R > but R xx = 0, for example, 
then R xx will first grow and then decay. 

4.3 Turbulent Lorentz forces and small-scale 
dynamo action 

If the fluid is initially at rest and a random magnetic field 
is introduced, the Lorentz forces will immediately drive tur- 
bulent motions. This is the origin of the term C3, which 
appears as a source for Rij and a sink for Mij. Note that 
a random magnetic field that has a dominant a;-component, 
say, drives motions predominantly in the :r-direction through 
the B ■ VB force. 7 For this reason the source for Rij has the 
tensorial form of My. 

The term C4 is a source for Mij and a sink for Rij . It can 
be considered as the effect of small-scale dynamo action. Any 
turbulent motion acts as a small-scale dynamo and causes 
a growth of turbulent magnetic energy, at least when the 
field is weak. This occurs through the B ■ Vu term in the 
induction equation, and therefore the source for My has the 
tensorial form of Rij. 

4.4 Energetics, equipartition and realizability 

The evolution of the turbulent kinetic and magnetic energy 
densities, ^R and |M, is governed by the traces of the equa- 
tions for Rij and My. Thus 

(d t + u t di)R = -2RijdjUi - C 1 L- 1 R Z/2 + CsL^M 3 ' 2 

-dL^R^M, (28) 

(d t + u t di)M = 2MijdjUi + C i L~ 1 R 1/2 M 

-(Ot + C^L^M 3 ' 2 . (29) 

7 This may seem paradoxical as the total Lorentz force is or- 
thogonal to B. However, only the solenoidal part of the Lorentz 
force drives motion in an incompressible fluid, or in a compress- 
ible fluid under Boussincsq conditions. The remaining gradient 
part is compensated by a pressure gradient. 



The total energy density satisfies 

(d t + Uidi)(±R + \M) = (M^ - Rij)djUi 

-C^L-^R 3 ' 2 - CsZT^M 372 . (30) 

The first term on the right-hand side represents the extrac- 
tion of shear energy from the mean flow by the total tur- 
bulent stress Mij — Rij. The other two terms, Ci and C5, 
represent the turbulent dissipation of kinetic and magnetic 
energy, and can be thought of as operating through the for- 
mation of vortex sheets and current sheets, respectively. 

Terms C3 and C4 transfer energy between kinetic and 
magnetic components. The net rate of transfer from kinetic 
to magnetic energy is 

(C 4 R 1/2 - Co,M xl2 )L- x \M (31) 

and is positive or negative according to whether the kinetic 
or magnetic energy dominates. There is therefore a tendency 
towards 'equipartition' in the ratio M/R — (C4/C3) 2 . 

An important requirement of the model is that the 
kinetic and magnetic energy densities should remain non- 
negative, not least because square roots appear in the model. 
In fact there is a stronger requirement, mentioned above, 
because the Reynolds and Maxwell tensors are by definition 
positive semi-definite tensors. In the Appendix it is shown 
that the positive semi-definite nature of Rij and Mij is pre- 
served by the model, guaranteeing that they are realizable 
by genuine velocity and magnetic fields. 



5 NON-LINEAR OUTCOME IN THE 
SHEARING SHEET 

In this section I examine the predictions of the model for 
the incompressible shearing sheet. The angular velocity and 
velocity gradient are fixed by the parameters of the shearing 
sheet. The Reynolds and Maxwell tensors may be assumed 
to be independent of position, corresponding to homoge- 
neous turbulence. The model then reduces to a non-linear 
dynamical system describing the purely temporal evolution 
of the stress tensors. The non-trivial fixed points of the 
dynamical system represent possible states of statistically 
steady turbulence. 

In numerical calculations I adopt the fiducial parame- 
ters Ci = C2 = C3 = C4 = Cb = 1, although other val- 
ues have been experimented with. The fact that this simple 
choice gives reasonable results indicates that no fine tuning 
is required for the model to behave physically. Ideally the 
parameters should be calibrated by comparison with suit- 
able numerical simulations. 

5.1 Freely decaying turbulence without rotation 
or shear 

With Q = A = there is nothing to initiate or sustain 
turbulence. Any turbulence introduced into the system de- 
cays, tending to isotropy as is does so. For freely decaying 
isotropic turbulence (Rij = ^RSij, etc.), the model states 
that 

R = -C 1 L~ 1 R 3/2 + CzL~ 1 M a/2 -C 4 L~ 1 R 1/2 M, 
M = C 4 L- 1 R 1/2 M - (C S + C 5 )L- 1 M 3/2 , (32) 
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Figure 1. Phase portrait for freely decaying isotropic turbulence 
in the absence of rotation or shear. Fiducial parameters Ci = 
C3 = Ci = C5 = 1 are adopted. The diagram is self-similar 
and therefore the scale is arbitrary. All trajectories tend towards 
the only fixed point, R = M = 0. During the decay there is a 
tendency towards cquipartition. 



Figure 2. Evolution of the Reynolds tensor in plane Couettc 
flow, starting from a small hydrodynamic perturbation without 
a magnetic field. Fiducial parameters C\ = C2 = 1 are adopted. 
The solution tends towards the stable fixed point representing a 
state of statistically steady hydrodynamic turbulence. The com- 
ponents R xz and R yz are decoupled from the others, and tend to 
zero. 



where the dot denotes differentiation with respect to time. 
A typical phase portrait is shown in Fig. 1. All trajectories 
tend towards the only fixed point, R = M = 0, but do not 
reach it in a finite time. During the decay there is a tendency 
towards equipartition. 



5.2 Non-rotating shear flow without a magnetic 
field 

The next simplest case to be examined is that of a non- 
rotating shear flow (plane Couette flow, O = 0) with Mij — 
throughout. This case has been studied in the laboratory 
(with rigid boundaries) and in numerical simulations with 
the periodic boundary conditions of the shearing box (Pumir 
1996). At high Reynolds number, sustained hydrodynamic 
turbulence develops. The closure model ought to describe 
this behaviour. 

As described in Section 3.3, in the absence of rotation 
the system is Rayleigh-neutral and the linear dynamics of 
Rij allows for algebraic growth. The outcome is determined 
by the non-linear terms. The system has only two fixed 
points, one of which is the algebraically unstable trivial so- 
lution Rij — 0. The non-trivial fixed point represents a state 
of steady hydrodynamic turbulence and is given by 



Rx 

Ry 

Rxv = 



R 



( — — — \ 

VC1 + C2/ 

( 3C 1+ C 2 \ 1 
VC1+C2 J 3 ' 



R, 



Ci 
ALA 

Ryz ' 



R 



3/2 



(33) 



R = 



C 2 



Ci(Ci+C 2 ) 2 



(34) 



The existence of this non-trivial solution depends on both 
the shear and the return-to-isotropy coefficient C2. The fixed 
point is stable for all parameter values and is universally 
attracting. 

A typical result of numerical integration of the time- 
dependent equations is shown in Fig. 2. The integration is 
started close to the unstable fixed point Rij — by intro- 
ducing a small positive value of R xx only. Initially algebraic 
growth leads to a rapid approach to the stable fixed point. 



5.3 Non-rotating shear flow with a magnetic field 

The dynamics changes when a magnetic field is introduced, 
because a turbulent flow acts as a small-scale dynamo. The 
non-trivial fixed point representing a state of steady hydro- 
dynamic turbulence is unstable to a magnetic perturbation. 
A new, stable fixed point appears, corresponding to a state 
of steady MHD turbulence. 

A typical result of numerical integration of the time- 
dependent equations is shown in Fig. 3. The integration is 



started close to the trivial solution Ri 



Mi. 



by in- 



with 



troducing small positive values of R xx and M xx only. The 
system approaches the stable fixed point representing mag- 
netized turbulence. 



5.4 Rotating shear flow without a magnetic field 

If one seeks a non-trivial fixed point representing a state of 
steady hydrodynamic turbulence in rotating plane Couette 
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Figure 3. Evolution of the Reynolds tensor (left) and Maxwell tensor (right) in plane Couctte flow, starting from small hydrodynamic 
and magnetic perturbations. Fiducial parameters C'i = C'2 = C'3 = C'4 = C5 = 1 are adopted. The solution tends towards the stable 
fixed point representing a state of statistically steady MHD turbulence. The solution approached in Fig. 2 is no longer attracting. The 
components R xz , Ryz, M xz and M yz arc decoupled from the others, and tend to zero. 



flow, one obtains the formal solution 
, 3Ro' 1 C\ + C 2 \ 1 



Rxx 
Ryy : 

Rzz l 

Rxy 

Rxz : 



C1 + C2 I 3 
3(l-Ro- 1 )Ci + C 2 



( C 2 

VCi + 



d + C 2 



d 

ALA 

Ryz ' 



R 



C*2 

3/2 



with 



R 



-6Ro _1 (Ro _1 - l)Ci + C 2 



Ci(Ci +C2) 2 



-L 2 A 2 , 
3 



(35) 



(36) 



which is the generalization of equation (34) for the non- 
rotating case. This solution is only meaningful, however, if 
Rij is positive semi-definite, and this requires 

§^ > GRo-^Ro- 1 - 1), (37) 
Ci 

as illustrated in Fig. 4. (In the case of equality, the solution 
is trivial.) Where it exists, this solution appears to be stable. 

For the Keplerian case Ro = 3/4 relevant to a stan- 
dard accretion disc, steady hydrodynamic turbulence is pos- 
sible only when C2/C1 > 8/3. Although this is possible in 
principle, such a value appears improbable a priori. When 
C2/C1 = 1, for example, steady turbulence is possible for 
—0.145 < Ro -1 < 1.145. For rotation laws of the form 
oc r~ 9 , this requires q > 1.746. 

Numerical integration of the time-dependent system 
confirms that, when the system is perturbed from Rij — 0, 
it tends towards the stable solution where this exists. Other- 
wise it exhibits decaying epicyclic oscillations about Rij — 0. 



5.5 Rotating shear flow with a magnetic field 

The introduction of a magnetic perturbation allows steady 
MHD turbulence to develop in a rotating shear flow, even 
when Rayleigh's stability criterion is amply satisfied, as in 
the Keplerian case Ro = 3/4 (Fig. 5). This occurs through a 
non-linear magnetorotational instability (non-linear because 
there is no imposed magnetic flux) and the magnetic field is 
sustained through a non-linear dynamo process. The model 
therefore reproduces in broad terms the finding of numerical 
studies of the non-linear evolution of the magnetorotational 
instability (e.g. Hawley et al. 1995). 



6 THE GOVERNING EQUATIONS IN 
COMPRESSIBLE FLOWS 

The assumption of an incompressible fluid is useful in re- 
ducing the problem to a minimal, although still formidable, 
complexity. There are two reasons why the model may not be 
suitable for immediate application to a compressible fluid. 
First, the mean density may be non-uniform, especially if 
the mean velocity field has a non-zero divergence. Second, 
the turbulence may be essentially transonic, changing the 
character of the motions and leading to greatly enhanced 
dissipation through shock formation. The second case is be- 
yond the scope of this paper and is unlikely to be important 
for the magnetorotational instability, at least in the part of 
an accretion disc where most of the mass resides. The first 
case, however, is important in a number of applications and 
can be treated in a simple way. 

In a compressible fluid it is more appropriate to define 
the mean Reynolds and Maxwell tensors as 



Rij - 

and 



{pUiUj) 



(38) 



Magnetorotational turbulent stresses 9 




Figure 5. Evolution of the Reynolds tensor (left) and Maxwell tensor (right) in a rotating shear flow with Ro = 3/4, starting from 
small hydrodynamic and magnetic perturbations. Fiducial parameters C'i = C'2 = C3 = C4 = C5 = 1 are adopted. The solution tends 
towards the stable fixed point representing a state of statistically steady MHD turbulence. The components R xz , Ryz, M xz and M yz 
are decoupled from the others, and tend to zero. 



Mij — ( — ) , (39) 
\ A*o / 

which have the dimensions of stress and differ from the ear- 
lier definitions by a factor of the density. It is convenient not 
to include the magnetic pressure perturbation in Mij, but 
to regard the mean magnetic stress as Mij — ^MSij. 

The first issue to consider is whether a divergence of 
the mean velocity field, diUi, should affect the evolution 
of Rij or Mij. This velocity divergence does not appear 
in the equation of motion, but does appear in the equa- 
tion of mass conservation and the induction equation in the 
forms dtp = • ■ ■ — pdiUi and dtBi — ■ ■ ■ — BidjUj respec- 
tively. This motivates the addition of terms in the forms 
dtRij = • • ■ — RijdkUk and dtMij = • • • — 2MijdkUk respec- 
tively. 

The second issue concerns the definition of the verti- 
cal length-scale L in a stratified disc, and the effect of the 
stratification on the vertical profile of the stress. In the in- 
compressible system the turbulence is homogeneous and the 
stress independent of z. Numerical simulations of stratified, 
isothermal accretion discs (Miller & Stone 2000) suggest 
that the stress is also stratified, being roughly proportional 
to the density (or pressure), and in many applications it 
is analytically convenient to assume that the stress scales 
with the pressure. This can be achieved in a natural way 
by identifying the vertical length-scale as L = c s /Q. z , where 
c s = (p/p) 1 / 2 is the isothermal sound speed and Q, z the ver- 
tical oscillation frequency, which measures the curvature of 
the gravitational potential and is equal to Q in a Keplerian 
disc. In an isothermal disc L is then equal to the Gaussian 
scale- height H. 

Accordingly, the system of equations recommended for 
a compressible accretion flow is as follows, written in an in- 
ertial frame of reference. The equation of mass conservation, 



d t p + di(pui)=0. (40) 
The equation of motion, 

p(d t +Ujdj)ui = -pdi$ - di(p+\M)+dj(Mij - Rij). (41) 

The equation for the Reynolds tensor, 

(dt + Ukdk)Rij + Rikdkitj + RjkdkUi + RijdkUk = 

n zP - 1/2 [-CiR^Rij - C 2 R 1/2 (R tJ - \RS tJ ) 

+C :i M 1/2 M t:i - C A R~ ll2 MR,j\ . (42) 
The equation for the Maxwell tensor, 
(d t + u k d k )Mij - M ik dkUj - M jk dkUi + 2Mijd k u k = 

n z p- 1/2 [C4.R~ 1/2 MRij - (C 3 + C 5 )M 1/2 M tj ] . (43) 
The thermal energy equation, 

P T(d t +uA)s = n zP - 1/2 (iCii? 3/2 + iC 5 M 3/2 ) -diFi.(44) 

Here T is the temperature, s the specific entropy and Fi 
the radiative energy flux. The positive contributions on the 
right-hand side represent turbulent heating through viscous 
and resistive decay. 

The total energy is then exactly conserved in the form 

9 t [p(i« 2 + $ + e) + ii?+iM] 

+di [p(^u 2 + $ + e)ui+pui + (\R + M)m 

+(R lJ - Mij) U j + Fi] = 0, (45) 

provided that the gravitational potential <3> is indepen- 
dent of t. Here e is the specific internal energy, such that 
de = T ds — pd(p _1 ). The existence of this conservation law, 
with positive-definite turbulent energy densities and heating 
rates, implies a certain self-consistency in the equations of 
the model. The terms that were added in passing to the 
compressible model are required to have the form that they 
do in order that energy be conserved. 
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Figure 4. Stability diagram for rotating plane Couette flow in 
the absence of a magnetic field. The two parameters of the system 
are the inverse Rossby number Q/A and the ratio C2/C1, which 
measures the tendency to return to isotropy. When the system 
is Raylcigh-unstable (0 < Q/A < 1), or when the isotropizing 
tendency is sufficiently large, a stable state of steady hydrody- 
namic turbulence is permitted. Otherwise the system exhibits de- 
caying epicyclic oscillations when perturbed. The Keplerian case 
Q/A = 4/3 is shown with a dotted line. Steady turbulence is 
permitted only if C2/C1 > 8/3. 

It is worth remarking that when Ci = • • • = C5 = 
and Rij = 0, these equations reduce exactly to those of ideal 
MHD if one identifies Mij as the Maxwell tensor BiBj / po of 
an arbitrary large-scale magnetic field £?; advected by the 
mean flow u;. For the magnetic terms in the equation of 
motion are precisely the Lorentz force of such a magnetic 
field, and the equation for the Maxwell tensor is equivalent 
to the induction equation of ideal MHD. Furthermore, the 
magnetic energy density, = B /2po, and the Poynting 
flux, Mui - MijUj = (B 2 Ui - BiBjUj)/no = (E x B) t /p , 
can both be identified in the energy equation. 



7 APPLICATIONS AND IMPLICATIONS 
7.1 Character of the equations 

The governing equations set out in the previous section usu- 
ally have the formal character of a hyperbolic system, indi- 
cating that information is propagated in a causal manner at 
finite speeds. To see this, consider the behaviour of infinites- 
imal wavelike perturbations of an arbitrary solution of the 
equations. Let the perturbations have the form of a rapidly 
varying phase factor, &cp[vp(x, t)], multiplied by slowly vary- 
ing functions of x and t. The wavenumber k — Vip and fre- 
quency lo = — d t ip are assumed to be such that |fe| _1 and 
I Cij I 1 are small compared to the typical length-scales and 
time-scales of the basic state. Although the model is not 
necessarily valid in this regime, such an analysis serves to 
uncover the formal mathematical structure of the equations. 



By linearizing equations (40)-(44) and eliminating per- 
turbations in favour the velocity perturbation i4, one ob- 
tains the algebraic eigenvalue problem 

p(u> — kjUj) 2 u'i = (Rjk + Mjkjkjkku'i 
+ (pv'jkikj + 2/i',, A' In': 

+(Mkikj - M lk kjk k - Mj k kik k )u'j (46) 

for the local WKB modes of the system, where v 2 — 
(dp/dp) s is the square of the adiabatic sound speed. As in 
compressible MHD, the WKB dispersion relation has 6 fre- 
quency eigenvalues u> for any choice of the real wavevector k. 
If any of the eigenvalues has a non-zero imaginary part, the 
system experiences a local instability. Otherwise the eigen- 
values are all real and the group velocities dui/dki are in- 
dependent of |fe|, indicating that the wave propagation in 
this limit is anisotropic but non-dispersive. The system of 
equations is then formally hyperbolic. 

A full investigation of the dispersion relation is difficult 
but the following observations may be made. First, when 
Rij — Mij = 0, the dispersion relation is the standard one 
of compressible hydrodynamics, and information is propa- 
gated at the sound speed relative to the mean flow. Second, 
when Rij — but Mij / 0, the system can be shown to 
be hyperbolic. The dispersion relation is related to that of 
compressible MHD. Third, in an incompressible fluid, modes 
with kiu'i =fc are eliminated as their speed of propagation is 
infinite. The remaining part of the dispersion relation, given 
by the first line of equation (46), describes the propagation 
of transverse modes, similar to Alfven waves, at finite speed. 
Finally, there are circumstances in which the system fails to 
be hyperbolic because the dispersion relation indicates a lo- 
cal instability, but this appears to be atypical. 

7.2 Stratified shearing sheet 

In a stratified shearing sheet, magnetorotational turbulence 
would develop according to these equations just as in an 
incompressible shearing sheet, except that all stress compo- 
nents would be proportional to p/Q. 2 rather than L 2 . The 
total shear stress M xy — R xy can be compared with that 
corresponding to an effective viscosity p — ap/Q. For the 
fiducial parameters in a Keplerian disc, this gives a w 0.081. 
If all the parameters C; are scaled by a constant factor, this 
value of a scales as 1/Cf . 

7.3 Non-Keplerian rotation 

When the Rossby number Ro = A/fl of the shearing sheet 
is varied in the range < Ro < 1 in which the system is 
Rayleigh-stable but magnetorotationally unstable, the total 
shear stress M xy — R xy in the steady turbulent state is not 
simply proportional to the shear rate 2A, at fixed angular 
velocity. This is shown in Fig. 6, and demonstrates one as- 
pect in which the present model differs significantly from a 
viscous representation of the turbulence. 

7.4 Damping of the two-dimensional wave 

The local dispersion relation for axisymmetric density waves 
in a two-dimensional disc model is 
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Figure 6. The non-linear relation between the total shear stress 
in steady magnetorotational turbulence and the shear rate, at 
fixed angular velocity. Fiducial parameters C\ = C2 = C3 = 
C4 = C5 = 1 are adopted. 

ui 2 = k 2 + v 2 k 2 , (47) 

where k is the radial wavenumber (Goldreich & Tremaine 
1979). The same mode can be identified within the present 
model il one considers a three-dimensional compressible 
shearing sheet without vertical gravity and stratification. To 
avoid complications it is convenient to replace the thermal 
energy equation by the isothermal condition p — c 2 p, c s = 
constant. Equations (40)-(43), linearized around a basic 
state consisting of steady, homogeneous magnetorotational 
turbulence, then admit solutions proportional to exp(ifca; — 
iut), with no perturbations of (u z , R xz , R yz , M xz , M yz ). In 
the absence of turbulent stresses, the dispersion relation (47) 
is recovered. A numerical solution of the dispersion relation 
in the presence of turbulent stresses indicates that the two- 
dimensional wave is damped in a Keplerian disc when the 
fiducial parameters C» = 1 are adopted (Fig. 7). The real 
part of the frequency is somewhat greater than predicted by 
equation (47), and agrees much better if v B is replaced by the 
'magnetosonic speed' (c s + M/p) 1 ^ 2 . Interestingly, the imag- 
inary part corresponds quite closely with a viscous damping 
rate vk 2 if the effective viscosity is v = ac 2 /fl with a w 0.06. 
This is remarkable when it is considered that there are no 
viscous or diffusive terms in the governing equations. 

7.5 Damping of shearing epicyclic motions 

When a stratified disc is warped, horizontal pressure gradi- 
ents are introduced that drive epicyclic motions with an am- 
plitude proportional to the distance z above the mid-plane 
(Papaloizou & Pringle 1983). The amplitude and phase of 
these oscillations are critical to the propagation and damp- 
ing of the warp. As discussed by Torkelsson et al. (2000), 
the dynamics of these shearing epicyclic motions, and in 
particular their damping time-scale, can be studied within 
the shearing sheet. 



Consider equations (40)-(43), linearized around a basic 
state consisting of steady, homogeneous magnetorotational 
turbulence, in which the perturbations depend only on z and 
t. It is possible to arrange for all perturbations to vanish 
other than those of (u x ,u y , R xz , R yz , M xz , M yz ). Then 

p(d t u' x - 2Ctu' y ) = 8 Z {M' XZ ~R XZ ), (48) 

p [dtu' y + 2(0, - A)u' x ] = d z {M' yz - R' yz ), (49) 

dtR' xz — 2Q,R' yz + R zz d z u' x — 

n z p- 1/2 [-(Ci + C 2 )R 1/2 R XZ + C 3 M 1/2 M' XZ 
-C A R~ 1/2 MR' XZ ] , (50) 

d t R' yz + 2(0, - A)R' XZ + R zz d z u y = 

V zP - 1/2 [-(Ci + C 2 )R 1/2 R' yz + C?M li2 M' yz 
-C A R~ 1/2 MR' yz ] , (51) 

dtM' xz - M zz d z u x = 

n z p- 1/2 [C 4 R- 1/2 MR' XZ - (C 3 + C 5 )M 1/2 M XZ ] , (52) 

d t M' yz + 2AM' XZ - M zz d z u' y = 

n z p- 1/2 [C 4 R~ 1/2 MR yz - (C 3 + C 5 )M 1/2 M yz ] . (53) 

These equations admit solutions in which all perturbations 
are proportional to exp(— iu>t), while u' x and u' y are propor- 
tional to 2 and R' xz , R' yz M' xz and M' yz are proportional to 
p. The six eigenvalues ui, corresponding to three physically 
distinct modes, follow from an algebraic eigenvalue prob- 
lem that is easily solved numerically. Two of these modes 
are strongly damped, while the other represents a slowly 
damped shearing epicyclic motion. For a Keplerian disc and 
fiducial parameters d = 1, the damping rate is 0.0254^, 
the same that would be obtained from an alpha viscosity if 
a = 0.0254. This is smaller than the value of 0.081 required 
to account for the shear stress M xy — R xy . 

8 COMPARISON WITH PREVIOUS WORK 

8.1 Comparison with analytical models by other 
authors 

The historical development of closed model equations for the 
stress tensor in a turbulent fluid is described in the review 
article of Speziale (1991). Early studies were based on the 
concepts of eddy viscosity and mixing length introduced by 
Boussinesq and Prandtl. A systematic investigation of the 
equation for the mean Reynolds stress was initiated in the 
1950s. Out of this arose (among others) the K-e model, 
which is widely used in engineering applications despite its 
numerous deficiencies. More recent approaches have aimed 
at a greater fidelity to experimental and numerical results, 
typically by elaborating algebraic models of the pressure- 
strain correlation and trying to fix the coefficients therein. 

The present work is related to these Reynolds-stress clo- 
sure models in that I have worked with the exact equations 
for the mean Reynolds and Maxwell tensors, retaining the 
form of the linear terms and proposing closures for the non- 
linear terms. However, I have tried to take a fundamentally 
different approach to the non-linear effects by recognizing 
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Figure 7. Dispersion relation of the two-dimensional wave in a three-dimensional, compressible shearing sheet without stratification. 
Keplerian rotation and isothermal perturbations are assumed, and fiducial parameters C\ = C'2 = C3 = C'4 = C5 = f are adopted. Left: 
real part of the frequency (solid line), compared with the standard dispersion relation (47) (dashed line). Much better agreement is found 
when the magnctosonic speed is used in place of the sound speed (dotted line). Right: damping rate of the wave (solid line), compared 
with the expectation based on a viscous damping corresponding to a = 0.06 (dashed line). 



the essential role of the length-scale L in the saturation 
of the turbulence and by identifying each non-linear term 
with a known physical process. In the absence of a mag- 
netic stress, the equations adopted here for the Reynolds 
tensor are considerably simpler than those advocated by, 
e.g., Speziale (1991), yet they do ensure realizability, allow 
for the return to isotropy, and may well give a more ac- 
curate representation of the linear and non-linear stability 
properties of rotating shear flows. 

The application of Reynolds-stress closure models to ac- 
cretion discs has been proposed in a series of papers by Kato 
(1994), Kato & Inagaki (1994) and Kato & Yoshizawa (1993, 
1995, 1997). In some of these papers the Maxwell stress is 
also modelled. While these studies have a fair amount in 
common with the present approach, some important differ- 
ences must be emphasized. In the absence of a magnetic 
stress, Kato & Yoshizawa (1997) have argued that steady 
hydrodynamic turbulence may be sustained in a Keplerian 
accretion disc. In fact, their closure model is based on one 
by Launder, Reece & Rodi (1975), which predicts stabil- 
ity for a Keplerian shear flow (see Speziale 1991); only by 
modifying one of the terms were Kato & Yoshizawa able to 
find a steady turbulent state. More importantly, the work 
of Kato and co-workers does not address the issue of the 
non-linear saturation of hydrodynamic or MHD turbulence. 
Although it predicts the relative magnitudes of the various 
components of the Reynolds (and Maxwell) tensors, it does 
not predict the magnitude of the turbulent energy in the 
saturated state. 

Some further connections with the astrophysical liter- 
ature can be made. Schramkowski et al. (1996) devised a 
kinematic mean-field theory for the evolution of the Maxwell 
tensor in a turbulent accretion disc, as an alternative to 
studying the mean magnetic field itself. Their theory is co- 



variant and includes the interaction of the Maxwell tensor 
with the mean velocity field. In addition to the 'alpha' and 
'beta' effects of standard mean-field electrodynamics, they 
emphasized the role of the 'gamma' effect by which a tur- 
bulent flow amplifies small-scale magnetic energy. This can 
be related to the C4 term of the present model. 

Soon after Balbus & Hawley (1991) first wrote on the 
magnetorotational instability in accretion discs, Tout & 
Pringle (1992) presented a simple model of a magnetic dy- 
namo cycle that uses the differential rotation to generate 
azimuthal field from radial field, the Parker instability to 
convert azimuthal field into vertical field, and the magne- 
torotational instability to regenerate radial field from verti- 
cal field, while magnetic reconnection provides dissipation. 
In effect, they wrote down a non-linear dynamical system 
for the three field components, which could be thought of as 
either mean or RMS values. The present model is ultimately 
based on a similar idea, although it treats the Maxwell stress 
covariantly as a single tensor object. The Parker instability 
(magnetic buoyancy) does not play a role here, and the mag- 
netorotational instability is an outcome of the model rather 
than an ingredient. 

8.2 Comparison with a viscoelastic model 

In recent work on the dynamics of eccentric discs I ques- 
tioned the use of a viscous model of the turbulent stress in 
an accretion disc (Ogilvie 2001). The assumption that the 
stress is linearly and instantaneously related to the rate of 
strain typically leads to an instability whereby an initially 
circular disc develops eccentricity on short radial length- 
scales (see also Lyubarskij, Postnov & Prokhorov 1994); in 
fact, this is essentially the same effect as the viscous over- 
stability of axisymmetric waves discovered by Kato (1978). 
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In order to take account of the non-zero relaxation time r 
of the turbulence, I proposed the use of a viscoelastic model 
for the Maxwell stress tensor, 

Mij + rVMij = ^(diUj + djUi) + (/Ub - |/tt) (d k u k )Sij, (54) 

where fi is the effective shear viscosity, /i b the effective 
bulk viscosity, and T> is a time-derivative operator such that 
VMij is equal to the left-hand side of equation (43). 

This model is a compressible version of the upper- 
convected Maxwell fluid of non-Newtonian fluid dynamics, 
and is motivated by the fact that the Maxwell tensor sat- 
isfies the equation DMij — exactly in ideal MHD. The 
analogy between viscoelasticity and MHD has been devel- 
oped in some detail by Ogilvie & Proctor (2002), who show 
that the magnetorotational instability is also found in vis- 
coelastic Couette flow. The introduction of the relaxation 
term, with r being comparable to the dynamical time-scale 
f^ 1 , has a profound influence on the dynamics of eccentric 
discs and tends to remove the viscous overstability. 

The viscoelastic model was intended to introduce the ef- 
fect of the relaxation time in the simplest possible way, and 
generalizes the conventional viscous model by the introduc- 
tion of a single new parameter, the Weissenberg number fir. 
It is not a true magnetorotational model because it neglects 
the Reynolds stress and relies on an effective viscosity to 
generate a non-zero Maxwell stress. In contrast, the model 
proposed in the present paper does not require any artificial 
source terms because the stresses arise through an instabil- 
ity of the trivial state and through a cooperative non-linear 
interaction between the Reynolds and Maxwell stresses. The 
cost paid for this is the need for approximately twice as many 
terms, variables and parameters in the model. Nevertheless, 
the viscoelastic model can be compared closely with equa- 
tion (43) if one identifies r" 1 = (C 3 + C 5 )n z (M/p) 1/2 and 
understands that the Ca term, in acting as a source for My , 
is analogous to the effective viscosity. 

8.3 Comparison with numerical simulations 

The behaviour of the present model in the shearing sheet, as 
set out in the various parts of Section 5, agrees at least qual- 
itatively with the results of existing numerical simulations. 
For a non-rotating shear flow without a magnetic field, one 
may compare with the simulations by Pumir (1996), which 
demonstrate the development of statistically steady, homo- 
geneous and anisotropic turbulence. A detailed comparison 
with the mean stresses obtained in numerical simulations, 
preferably in the limit of horizontally extended shearing 
boxes L x ,L y ^> L z , would test the accuracy of the model 
and constrain the values of the parameters. 

Rotating shear flows without a magnetic field have been 
studied numerically by Hawley, Balbus & Winters (1999). 
They investigated rotation laws of the form Q oc r~ q , find- 
ing that Keplerian shear flows (q — 3/2) are hydrodynam- 
ically stable while constant-angular momentum shear flows 
(q — 2) are non- linearly unstable and become turbulent. 
They tried to locate the critical value of q at which hydro- 
dynamic turbulence could just be sustained, and found it to 
be close to 1.95, although this may depend to some extent 
on the resolution of the numerical method and the initial 
conditions adopted. These results are in qualitative agree- 
ment with Fig. 4, which also suggests that hydrodynamic 



turbulence can be sustained only in flows that are Rayleigh- 
unstable or only just Rayleigh-stable. 

The present model predicts the development of steady 
MHD turbulence in a Keplerian shear flow, much as seen in 
numerical simulations by Hawley et al. (1995), although it is 
more directly applicable to calculations without an imposed 
magnetic flux, such as those of Brandenburg et al. (1995) 
and Stone et al. (1996). The initial magnetic field in the 
latter simulations is well ordered, although of zero mean, 
and gives rise to an initial phase of exponential growth that 
is not captured in the model. Nevertheless, the properties 
of the saturated stresses appear to be in good qualitative 
agreement. A detailed numerical comparison would be in- 
appropriate, because unfortunately there is no consensus as 
to the 'correct' absolute or relative magnitudes of the mean 
stress components in a Keplerian shear flow. For example, 
Brandenburg et al. (1995) and Stone et al. (1996), although 
adopting identical initial conditions and numerical resolu- 
tions, obtain widely differing quantitative results; in partic- 
ular, they disagree on whether the turbulent kinetic energy 
or the turbulent magnetic energy is greater. It is unclear 
whether this discrepancy is attributable to the difference in 
boundary conditions or to a lack of numerical convergence. 

Abramowicz et al. (1996) have simulated the magne- 
torotational instability in non-Keplerian shear flows, having 
in mind the application to the inner parts of accretion discs 
around black holes. They found a non-linear relation be- 
tween the shear stress and the shear rate, at fixed angular 
velocity, in good agreement with Fig. 6. 

Interestingly, the present model also allows for sustained 
MHD turbulence in the case of negative Rossby number, 
A/Q < 0, which corresponds to a situation in which the 
angular velocity increases outwards. This regime may ap- 
pear astrophysically improbable but may be relevant to the 
boundary layer between a star and a disc. While bound- 
ary layers have been simulated in recent global studies (Ar- 
mitage 2002; Steinacker & Papaloizou 2002), the case of neg- 
ative Rossby number does not appear to have been explored 
in shearing-box simulations. There is no linear magnetoro- 
tational instability in this case (except when the Hall effect 
is important, cf. Balbus & Terquem 2001), but the present 
model suggests that a turbulent MHD state may neverthe- 
less exist and is accessed by a non-linear instability of the 
basic state. The energetic arguments presented by Balbus 
& Hawley (1998) appear to allow for this possibility, and 
it would be valuable to test this hypothesis in numerical 
simulations. 

In one study, Hawley, Gammie & Balbus (1996) re- 
peated one of their typical shearing-box simulation of mag- 
netorotational turbulence, but omitted the Coriolis force, 
thereby converting the calculation into a local model of a 
non-rotating shear flow. They found that the magnetic en- 
ergy was not sustained, although the turbulent kinetic en- 
ergy continued to grow until the calculation was halted. 
This differs from the prediction of the present model in Sec- 
tion 5.3, that turbulent plane Couette flow acts as a small- 
scale dynamo and allows the magnetic energy to grow to 
a substantial fraction of equipartition with the turbulent 
kinetic energy, which itself reaches a saturated value. In- 
deed, the model is constructed on the basis that any three- 
dimensional turbulent flow at sufficiently high Rm, by virtue 
of its stretching property, tends to amplify a weak, small- 



14 G. I. Ogilvie 



scale magnetic field introduced into it, at least until the 
Lorentz forces can modify the flow. In plane Couette flow, 
the shearing background provides a further means of ampli- 
fying the magnetic energy. The numerical result of Hawley 
et al. (1996), that the combination of turbulence and shear 
has the opposite characteristic, is therefore surprising. One 
possible explanation for the discrepancy between the numer- 
ical result and the prediction of the model is that Rm is too 
small in the simulation to capture the small-scale dynamo. 
Alternatively, it may indicate a deficiency in the model. Fur- 
ther numerical studies may resolve this question. 

Regrettably, there are no published numerical studies of 
the propagation and damping of the two-dimensional wave 
in turbulent discs, to compare with the predictions of Sec- 
tion 7.4. This problem is of some importance, especially in 
view of its application to tidally generated waves in binary 
stars and protoplanetary systems, and would make a further 
test of the present model. 

However, Torkelsson et al. (2000) have already studied 
numerically the interaction of magnetorotational turbulence 
with the shearing epicyclic motions found in warped discs 
(cf. Section 7.5). They found that the epicyclic motions are 
damped, but at a rate approximately one-half that expected 
from the action of an isotropic viscosity (the viscosity being 
measured from the total xy-stvess). This is similar to the 
predictions of the present model as set out in Section 7.5. 



9 CONCLUSION 

The representation of the turbulent stresses in an accretion 
disc by the classical alpha viscosity model of Shakura & 
Sunyaev (1973) has been extremely successful in allowing 
the theory of accretion discs to develop even in the ab- 
sence of a proper understanding of the turbulence, and in 
making possible a quantitative comparison between theory 
and observations. However, the magnetorotational instabil- 
ity is now widely accepted to be the origin of the turbulence 
in sufficiently ionized discs (e.g. Balbus & Hawley 1998), 
and the alpha viscosity model fails to describe numerous 
aspects of this process. In this paper I have introduced a 
new analytical model that aims to represent more faith- 
fully the dynamics of magnetorotational turbulent stresses 
and bridge the gap between analytical studies and numer- 
ical simulations. Covariant evolutionary equations for the 
mean Reynolds and Maxwell tensors are presented, which 
correctly include the linear interaction with the mean flow. 
Non-linear and dissipative effects, in the absence of an im- 
posed magnetic flux and in the limit of large Reynolds num- 
ber and magnetic Reynolds number, are modelled through 
five non-linear terms that represent known physical pro- 
cesses and are strongly constrained by symmetry properties 
and dimensional considerations. 

The cooperation of these linear and non-linear effects 
explains, without any fine tuning, the development of sta- 
tistically steady, anisotropic turbulent stresses in the shear- 
ing sheet, a local representation of a differentially rotating 
disc, in agreement with numerical simulations. It reproduces 
other results seen in numerical studies: that purely hydrody- 
namic turbulence is not sustained in a flow that adequately 
satisfies Rayleigh's stability criterion; that the shear stress in 
magnetorotational turbulence in non-Keplerian discs is non- 



linearly related to the shear rate, at fixed angular velocity; 
and that the shearing epicyclic motions found in warped 
discs are damped by the turbulence at a slower rate than 
predicted on the basis of an isotropic viscosity. It also makes 
predictions for future simulations, such as the decay of the 
two-dimensional wave and the possibility of sustained tur- 
bulence in situations where the angular velocity increases 
outwards. The model has good formal properties, being typ- 
ically hyperbolic and therefore 'causal', guaranteeing the re- 
alizability of the stress tensors, and accounting satisfactorily 
for energy conservation. 

Only the linear part of the model is derived directly 
from the MHD equations, and therefore the predictions of 
the model cannot be regarded as rigorous deductions. In par- 
ticular, the non-linear hydrodynamic stability of circular Ke- 
plerian flow remains unproven. Nevertheless, even without 
an accurate calibration of its parameters, the model almost 
certainly performs more faithfully than the alpha viscosity 
model in a variety of circumstances, and therefore represents 
an advance over previous work. 

The general idea of seeking a model that connects the 
mean turbulent stress to changes in the mean velocity gra- 
dient can be related to one of the classical problems of con- 
tinuum mechanics, which is to find the constitutive rela- 
tion that characterizes a particular substance and allows 
the stress in a body to be related to its deformation history. 
From this perspective, the present model is reminiscent of 
certain classes of viscoelastic fluids. In fact, the coexistence 
of the Reynolds and Maxwell stresses, which are influenced 
by the mean flow in different ways, suggests a kind of com- 
posite viscoelastic material. The reason that a constitutive 
relation might be thought to exist at all for magnetorota- 
tional turbulence is that the physics of the magnetorota- 
tional instability is local, suggesting that the mean stress in 
a volume comparable to H 3 should indeed depend only on 
the recent deformation history of that parcel of fluid. 

In recent years, numerical simulations of accretion flows 
subject to the magnetorotational instability have been taken 
beyond the confines of the shearing box to study the effects 
of cylindrical geometry, the flow across the marginally sta- 
ble orbit around a black hole, or the boundary layer between 
a star and a disc (e.g. Armitage 1998; Hawley 2000; Haw- 
ley, Balbus & Stone 2001; Armitage 2002). These advances 
are to be welcomed, and represent significant computational 
achievements. However, I believe there is still, and will con- 
tinue to exist, a need for a description such as the one de- 
veloped in this paper. This is not just because global studies 
of realistic thin discs are well beyond the current capabili- 
ties of direct numerical simulations, but also because such 
a model will allow a variety of methods, analytical and less 
direct numerical ones, to continue to be used in studying 
accretion disc dynamics while taking seriously the nature of 
magnetorotational turbulence. It should be particularly use- 
ful in understanding the dynamics of warped, eccentric and 
tidally distorted discs, non-Keplerian accretion flows close 
to black holes, and a variety of time-dependent accretion 
phenomena. 
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APPENDIX A: REALIZ ABILITY 

The realizability constraint requires that the three eigen- 
values of Rij and the three eigenvalues of Mij remain non- 
negative under the evolutionary model. A full investigation 
of this issue would require a consideration of a large number 
of degenerate cases. For simplicity, I demonstrate here only 
the weaker result that no one of the eigenvalues can become 
negative while the others remain positive. 

It is therefore assumed that, in the initial condition, 
Rij and Mij are positive definite tensors. If either Rij or 
Mij were to develop a negative eigenvalue in the subsequent 
evolution, when the first such eigenvalue passed through zero 
the traces R and M (which are the sums of the eigenvalues) 
would still be positive. Therefore one may take R, M > 
below, without loss of generality. 

Define the quadratic forms P = RijXiXj and Q — 
MijYiYj, where X, and Y; are differentiable vector fields 
advected according to the time-reversible equations 

(d t + Ujdj)Xi - XjdiUj + 2e ikl Q, k X l = 0, (Al) 

(d t + Ujdj)Yi + Y 3 d,u 3 = 0. (A2) 

Then 

^ = -(Ci + C7 2 )i- 1 i? 1/2 P + ICzL-'R^XiXi 

+C a L" 1 M 1/2 M iJ X i Xj - C 4 L~ 1 R- 1/2 MP, (A3) 

2J2 = CiL^R^^MRijYiYj - (C* 3 + C^L' 1 M 1/2 Q,(M) 

where D/Dt = d t + Uidi is the Lagrangian time-derivative. 

In the initial condition, P > everywhere for all Xi 
and Q > everywhere for all Y;. According to equation 
(A3), DP/Di > whenever P = and Mij is positive def- 
inite, and therefore P cannot become negative. According 
to equation (A4), DQ/Dt > whenever Q — and Rij is 
positive definite, and therefore Q cannot become negative. 



